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(revision of January, 2000) and implemented in Fortran 77. The algorithm handles both c.m.s. and 
hadron collision kinematics. 



Introduction 1 

Jet finding algorithms are a central data-processing tool in 
high-energy physics. However, the definition of jets suffered 
from misinterpreted ambiguities, which resulted in several 
existing jet finding schemes (cone, recombination, . . .), each in 
several variations (for a review see e.g. [1]). This resulted in 
practical difficulties since each experiment tends to use its own 
jet definition, making it difficult to compare physical results. 

A systematic theoretical analysis of jet-related measure- 
ments from first principles was performed in [2], [3] where the 
central importance of the following two general requirements 
was pointed out: 

(a) stability of data processing algorithms with respect to 
small effects such as errors in experimental data, etc.; 

(b) compatibility with quantum field theory (the connection 
of generalized shape observables with the fundamental energy- 
momentum tensor was established in [4]). 

The theory of [2], [3] criticized the conventional view on 
jet-finding algorithms as inversion of hadronization and, 
instead, offered to regard jets as a tool of approximate descrip- 
tion of hadronic events — in full compliance with the first 
principles of physical measurements [2], [3]. This opened way 
for deriving a jet finding criterion entirely from first principles 
[5]. The criterion is derived from studying how physical infor- 
mation (represented by the mentioned fundamental shape 
observables) is distorted in the transition from particles to jets, 
and requiring that such a distortion is minimal. 

The criterion of [5] can be regarded as a cone algorithm 
rewritten in terms of thrust-like shape observables, which 
makes it potentially very attractive. 

In the practical aspect, however, the algorithmic usefulness 
of the definition [5] has remained somewhat problematic 
because it involved a minimization of a function within a 
bounded domain in the space of iVj ets X A^ particles dimensions 

with an additional restriction. With A^ partic i es = 100-^400 (the 



numbers typical for LHC) and iVj ets = 1 -r 6 , one deals with a 

huge dimensionality (up to 2000 and more), which makes the 
resulting optimization problem non-trivial. 

Fortunately, the analytical regularity of the criterion of [5] 
means there is a lot of information about the problem — 
enough to allow a quite efficient implementation of the corres- 
ponding jet-finding algorithm, apparently on a par with con- 
ventional algorithms. 

The purpose of the present Letter is to describe such an 
algorithm with emphasis on analytical, programming-language- 
independent aspects. The currently available version is written 
in Fortran 77 (adapted from a code developed in Component 
Pascal [7]), has been compiled on a number of platforms, and 
works for both c.m.s. (spherically symmetric) kinematics, and 
for hadron-hadron collisions (cylindrical kinematics). The 
algorithm reliably solves the minimization problem of [5] (i.e. 
finds the optimal configuration of jets) for a typical event in a 
small fraction of a second on a modest workstation. 

The source together with a practical description of interfa- 
ces and examples is available on the Web. We aimed at de- 
signing a generic robust algorithm and well-structured code 
with a sufficient number of interface hooks to allow one to 
perform further modifications in case of need. Here we just 
mention that the central subroutine is Q_minimize; it de- 
scends into a (local) minimum from a given initial configura- 
tion, and it can be used in different ways, depending on a spec- 
ific application, such as finding all local minima, etc.; cf. the 
discussion of construction of jet-based observables in [5]. 
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The criterion 2 

In the summary description given below we do not attempt 
to explain the meaning of, and motivations behind all the 
details of the definition; the reader is referred to [5]. 

The event P is represented as a collection of its "particles", 



For each jet one also defines a light-like Lorentz 4-vector 



P ={ E «>Pa} a 



2.1 



each particle described by its energy E a and direction p a , the 

interpretation of which two quantities depends on the 
kinematics: 

For spherical ( e + e~ ) kinematics, E a is the energy of the 
particle, and p a is the corresponding unit 3-vector (which can 
be represented e.g. by the two standard angles (p a ,8 a )• 

For cylindrical (hadron collisions) kinematics, E a is the 
transverse energy (also denoted as E la = E" ne sin0 a ), 
whereas p a is represented by the (pseudo) rapidity 
7] a = lncot(0 a 12), - oo < 7j fl < +oo (the polar angle 9 a is 
measured from the beam axis), and the transverse direction 
Pla ( umt 2-vector which is normal to the beam axis; it can be 
represented by (p a ). 

In either case, the energies E a are normalized so that 

E E„ = 1 . 2.2 

t—ia a 

With each pair E a ,p a , one associates a unique light-like 

Lorentz 4-vector p a ,p\ = . 

A configuration of jets Q is described in a similar fashion as 
a collection of jets' energies and directions: 



jets 
= 1 



2.3 



In the conventional algorithms, particles and jets are linked 
by indicating which particle belongs to which jet. The scheme 
of [5] allows more flexibility: one introduces the so-called 
recombination matrix z a j which describes the fraction of a-th 

particle which went into formation of the 7-th jet. z aj can be 

any real number < z a j < 1 such that 



Ljfz^l for any a. 



2.4 



The inequality corresponds to the fact that part of the particle's 
energy is allowed to not participate in jet formation. So it is 
convenient to introduce the quantity 



z a - x ~lLj z aj foran y a 



2.5 



Note that the conventional scheme corresponds to restricting 
z a j to the value 1 if the a-th particle belongs to the j-th jet, 

and to in the opposite case. 

The recombination matrix 7 . is the fundamental unknown 

UJ 

in our scheme. In particular, the quantities interpreted as jets' 
physical momenta p . are expressed as follows: 



1j' 1j 



called the jet's 4-direction , as follows: 



Spherical kinematics: 
qj=(lqj), 2.7 

where = pj / \pj\ is a unit 3-vector and pj is the space- 
like component of p . . 

Cylindrical kinematics: 
q - = (cosh 77 ., sinhr/ ; ., qj) , 2.8 

where qj- = pj/\pj\ is a unit 2-vector with pj- being the 
transverse component of pj , and 77 • determined from the 
relation 

Further, one defines two functions: 
Y[P,Q] = 2^.p..^., E soft [P,Q]=£ a Z a £ a . 2.10 

The latter is interpreted as the "soft energy" which does not 
take part in jet formation. 

The criterion is as follows. One chooses R > and defines 



£2 R [P, Q] = R 2 Y[P, Q] + E soft [P, Q] 



2.11 



Then one chooses a (small; say, 0.01) number o cut > and 



finds z a j which minimizes £2 and satisfies the restriction 



Sl R [P,Q] <co c 



2.12 



with a minimal number of jets. 

It turns out that this jet finding criterion is similar to 
conventional cone algorithms although the expression 2.11 is a 
shape observable that generalizes the thrust to any number of 
thrust (semi-) axes (see [5] for a detailed discussion). 
Correspondingly, the parameter R is similar to the cone radius 
of the conventional cone algorithm (the standard value 0.7 
often adopted in the conventional algorithms roughly corres- 
ponds to R = 1 in our case). However, jet shapes in our case 
are determined dynamically taking into account the global 
shape of energy flow of the event. 

The physical meaning of R is the maximal jet radius as 
measured by infinitesimally soft particles (i.e. such a particle is 
relegated to soft energy if it is farther than R from any jet's 
axis). 

The parameter co cut is analogous to the jet resolution 
parameter of conventional recombination algorithms — and, 
simultaneously, to the so-called f-cuts in conventional 
algorithms [8] because Eq.2.12 imposes an upper bound on 
soft energy. 

See [5] for a detailed discussion of all this. 



Pj = L a z aj p a . 



2.6 



D.Yu.Grigoriev and F.V.Tkachov, hep-ph/9912415 [corrected 2001-08-15] Implementation of optimal jet definition 



The algorithm 

It is convenient to treat the "soft energy" formally as a 
special 0-th jet and denote z Q j = z a ■ Then 



E^o% = 1 foran y «• 



3.1 



The AT^-dimensional region described by the restriction 3.1 is 



the standard N. -dimensional simplex. Thus the configuration 



jets 

z a j one has to find is a point in a iVj ets x iV part -dimensional 

region which is a direct product of iV part iVj ets -dimensional 
standard simplices. 

The algorithm we found to work well is a hybrid of the 
gradient method and a coordinate-by-coordinate optimization 
as well as a heuristic based on the experimental finding 
(a posteriori supported by some theoretical arguments) that the 
minimum tends to be located on configurations with the matrix 
elements z a j taking the values of either zero or one, which 
corresponds to vertices of the simplices 3.1. The algorithm 
consists in iteratively performing minimization with respect to 



-a] 



within the simplex 3.1 for each particle a. 



Gradient minimization within the standard simplex 3.2 
For simplicity denote n = iVj ets and let Zj = z a j ■ The 

corresponding n -dimensional vector is denoted as z. It is 
convenient to work in terms of finding a maximum rather than 
minimum. So we are going to find a maximum of the function 
F(z) = -Q. R (z) forz within a domain D described by simple 
linear inequalities ("the standard n -dimensional simplex"): 



Zj > 0, 



, < 1 ■ 



3.3 



Sums over j such as in 3.3 run over j = 1 n unless 

explicitly restricted. 

The simplest algorithm of maximum search is to start from 
a candidate point z and to find the next candidate point in the 
form 



zd 



3.4 



where z > is a number and d describes a direction in which 
F increases. We are not interested here in z (see however 
Sec. 4) and focus on finding d up to an overall scalar factor 
from the local properties of F (its first derivatives at z) taking 
into account that z+zd must remain within the boundaries of 

the domain D . 

First suppose z is an internal point of D . The function is 
locally represented as 



F(z + zd)=F(z) + z{f-d) + 0(d 2 ) 
where fj =dF(z)/dzj and 

(fd) = Z J f i d ] - 

A natural desire is to find the 
direction d in which the linear 
function (fd) increases fastest. 

However, in order to quantify such 
a desire, one has to define the 
notion of distance along each 



3.5 



3.6 




(/ d) = const 



d l =\ 



direction, and this involves an arbitrariness. For example, one 
could use euclidean distance and then the direction of fastest 
increase corresponds to a uniquely defined point of the unit 
sphere — but the choice of the euclidean metrics itself is not 
unique (cf. the figure). In general, any vector satisfying 
(/ d ) > can be made the direction of fastest increase for 

some euclidean metrics. The only general heuristic is to choose 
d in some simple way which is natural in the context of a 
specific problem. If the only information available is Eq.3.5 
then the usual choice is 



d=f. 



3.7 



We adopt this choice for internal points of D, whereas the 
mentioned arbitrariness is made use of in determining d in the 
case when z is on the boundary of D. 

Next suppose z is a point of the boundary of D for which 

Zj=0, jeB°; zj>0 j£B°; LjZj<l, 3.8 

where B° is a subset of {1, ... , n} . Then d must obey the 
following restrictions: 

dj > 0, j £ B° . 3.9 
A simple choice for such d is this: 

IF j e B° THEN dj := MAX(0, /•) ELSE dj := fj END . 3.10 

Consider the case when z sits on the front face of the 
simplex 3.3, i.e. 



3.11 



Then instead of 3.9 one has 

y.d,<0. 3.12 

Choose 1 < J < n and change coordinates by replacing dj 
with the new independent coordinate d Q : 

d o = -Lj d j' d i ■ Y. <>■■■■■ 

Note the useful symmetry between all the components 

d i , j = 0, ... ,n which is best seen from the relation 

J 



3.13 



3.14 



3.15 



In terms of 3.13, the restriction 3.12 takes the standard form 
d > . Re-express (f-d) in terms of the new independent 

coordinates (dj^j, rf ) : 
where 

fj*j=fj-fj, /o ="/./■ 3.16 
In terms of the coordinates (dj^j , rf ) , a valid direction of 
increase for (f d) is given by the following analog of 3.10: 



IF j = THEN dj = d Q :=MAX(0, / ) 

ELSIF j*J THEN d : \= f : END. 
J J J J 

Lastly, d f is found from 3.13. 



3.17 
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We can now consider the points of the boundary of D which 
satisfy the most general set of restrictions: 

Zj = 0, jeB°; Zj >0, j £ B° ; ZjZj<l. 3.18 
This can be represented in a symmetric fashion by introducing 

Zo = l-E,.z r 3.19 

Then let B = B° u {0} c {0, ... , n] , and Eq.3.18 becomes 
Zj=0, jeB; Zj >0, j£B. 3.20 

Note that B * {0 n] (no point on the boundary of the 

simplex can belong to all its faces simultaneously). Therefore, 
one can always choose J <£ B and consider d • , j * J as 

independent variables (here j can take the value 0). Then from 
3.15 one deduces the following prescription for choosing d: 

IF jeB THEN dj :=MAX(0, /•) 

ELSIF j*J THEN d , := f '■ END. 3.21 

J J J J 

Lastly, dj is found from 3.13. 

The above choices of d are not unique as seen from the 
arbitrariness in the choice of J. 



Formulas for derivatives 



3.22 



Here are the relevant formulas for derivatives of Y with 
respect to z a f 



-Y(M) = 2 



= 2p a qj+2 Pj 



\ a ' J 



3.23 



Only one term in the sum over j survives because the terms 
which correspond to j'* j depend only on z a y. 

To evaluate the second term one has to use specific 
expressions for each of the two standard kinematics. For 
spherical kinematics we obtain 

= -\pj\- 1 [(p r P a )-(p r q j )(q rPa )]^0. 3. 
For cylindrical kinematics we obtain: 



"3V 



24 



dz * \~ p j 



9z„ 



dz aj 



[ p°j sinh?7 ; . - pj coshr/j ] , 



where 



= El)E ±a {r la -r lj ). 



3.25 



3.26 



The resulting formulas are sufficiently simple not to present 
calculational difficulties in either case. 



Implementation 4 

We limit our discussion here to language-independent 
aspects. Specific interfaces and code examples are provided 
separately [6]. We only note that the design of algorithm 
(which required experimentation with data structures and 
interactive experimentation with the control parameters of the 
algorithm) was performed using the freely available RAD tool 
BlackBox Component Builder [7] based on an extremely 
simple, type-safe, object-oriented and efficient compiled 
language Component Pascal (of the distinguished 
Pascal/Modula-2/Oberon-2 pedigree [9]). The final algorithm 
turned out to be simple enough to allow a port to Fortran with 
some improvements resulting from experimentation with 
realistic test samples of events. It should be emphasized that 
the design of the algorithm would have been much harder 
without all the safety features and the stunning combination of 
power and simplicity of Component Pascal, and without the 
simplicity, high interactivity and GUI features of the BBCB. 

Concerning our test samples of events, we used the total of 
500 events generated using Jetset [10] for typical processes 
studied at CERN and FNAL. It was not our aim to arrive at any 
physical conclusions, and in fact the specific nature of events 
played practically no role because our algorithm is fairly 
generic, and its overall behavior is essentially insensitive to 
details of structure of events. The tests were performed only for 
numerical debugging, not any studies of physics. A final 
adjustment of some numerical parameters was made possible 
by a large-scale test run on a realistic event sample performed 
by Pablo Achard [11]. 

No comparison with conventional algorithms has been 
attempted yet (the situation has changed by 2001 — FT). 

The minimization scheme described above is easily and 
straightforwardly implemented using only static data structures 
(easily mapped to Fortran arrays), among which the central are 
the 2-dimensional array z a j and the 1 -dimensional array 

corresponding to the direction d. The total data size is 
determined by the size of z a j. The number of particles cannot 
exceed a few thousands, and the number of jets, a dozen or so. 
So the size of z aj is 0(1000) x O(10) x (8 Bytes) ~ O(W0K) . 

If each subarray z a j for fixed a is contained in a contiguous 
memory block then the internal loop (which corresponds to 
minimization with respect to one particle's parameters) always 
deals with 0(lK ) of contiguous data, which ensures a very 
good localization of the algorithm and therefore a fast 
performance. 

Concerning the ambiguity in the choice of d according to 
the formulas given in Sec. 3.2, we found it advantageous to 
perform maximization of the length of d (measured according 
to the simple norm \d \ = max 7=0 n | dj\ which is natural in 

the context of simplicial geometry) with respect to J (which is 
a free parameter in the above formulas). Such an optimization 
involves a small amount of well-localized data and code 
involving only very simple operations, resulting in a fast 
execution, whereas the resulting overall speedup proved to be 
significant. 

The choice of step length % (cf. 3.4) is determined by the 
experimental finding that the minimum tends to be located at 
the boundary of the simplex. So we find T from the require- 
ment that the new candidate minimum z + td for given z and 
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d should be located at the boundary of the simplex, and if this 
results in an increase in the value of £1, then % is iteratively 
divided by a constant factor ~3 until a minimum is found. 

An important technical implementation detail (the so-called 
"snapping") concerns how one deals with the boundaries of the 
simplex: if some z a j is small enough (i.e. z is close enough to a 
boundary of the simplex) then it is set to zero ("snapped" to 
the boundary). A similar snapping is used for the direction d. 
Such snappings are necessary because one needs to detect the 
situations when z is at the boundary and the direction runs 
parallel to the boundary. 

There are no difficulties with the termination condition: 
since the resulting minimum is located at the boundary of the 
region (we have not seen exceptions so far, and some analytical 
arguments seem to indicate that such exceptions can never 
occur), the minimum tends to "snap" to the boundary quite 
fast. Also, most particles find their jet pretty fast, and later 
iterations involve decreasingly smaller numbers of particles. 

The fact that the trajectory of the search tends to travel 
along the boundary of the region makes the algorithm similar 
to the well-known simplex algorithm of combinatorial 
optimization but in our case the algorithm is sped up by 
reliance on explicit analytical formulas to determine the 
direction of the next candidate minimum. 

Typically, the algorithm arrives at a (local) minimum from a 
randomly generated starting point in 0(100) iterations, in a 
fraction of a second on a Pentium II. 

In general, for a given event the criterion may have more 
than one local minimum. This is discussed in detail in [5]. 
To find the global minimum it proved sufficient to repeat the 
search a few times starting from new randomly generated con- 
figurations z a j- "A few times" depends on the character of 
one's events: 2 (or even 1) may be enough for most situations 
with hard jets, and 10 seemed to be sufficient for events corres- 
ponding to the LEP2 process e + e~ — > W^W" — > 4 jets. The 
number of repeated searches anticorrelates with the fraction of 
events for which the algorithm fails to correctly identify global 
minimum. This number is therefore tied to the overall preci- 
sion of the physical problem. The implementation also 
provides for an explicit specification of the initial configuration 
to allow e.g. output from a conventional algorithm to be used 
for that purpose. 

Further optimizations are possible by way of adding more 
intellect/memory to the algorithm (e.g. giving priority in 
minimization to some particles, or using special heuristics to 
reduce the number of repeated searches in situations where 



several local minima may occur) but we felt that in view of a 
good speed of the minimum search the additional complexity is 
not warranted at this stage. So we limited the design to a 
generic algorithm while providing the modularity to allow one 
to build such improved algorithms in case of real need. 

To summarize, the described implementation proves practi- 
cal feasibility of the jet definition of [5], and the developed 
software allows easy modifications to accommodate further 
data processing options described in [5]. Such options, while 
potentially important in specific applications from physical 
viewpoint, are not expected to require major changes of the 
described minimization algorithm. 
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